В качестве сырых данных возьмем анонимизированные данные о поездках желтых такси в Нью-Йорке.
Сырые данные о поездках в жёлтом такси можно найти на сайте TLC: https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page
Эти данные разбиты на файлы по месяцам.
Будем использовать информацию о поездках до мая 2016 года включительно, т.к они имеют более удобный вид для обработки.
Задача проекта — научиться предсказывать количество поездок на ближайшие сутки (24 часа) в каждом районе Нью-Йорка;
Для простоты мы поделим область Нью-Йорка сеткой 50 на 50. Поскольку в проекте нужно предсказывать 1 временной ряд, то мы выберем одну из центральных областей (под номером 1231), и будем предсказывать количество поездок именно для нее. Но в перспективе нужно работать со всеми ячейками.
Для того, чтобы это сделать, сырые данные необходимо агрегировать по часам и районам.
Агрегированные данные будут представлять собой почасовые временные ряды с количеством поездок из каждого района.
За основу проекта брался мой недоделанный проект для Courser'ы, поэтому возможно у кого то будет такая же идея.
import pandas as pd
import numpy as np
from scipy import stats
data = pd.read_csv("yellow_taxi_data\yellow_tripdata_2016-05.csv" ,sep=',', header=0,
parse_dates=['tpep_pickup_datetime', 'tpep_dropoff_datetime'], infer_datetime_format=True,
dayfirst=True, engine='c', encoding='cp1251')
data.head()
| tpep_pickup_datetime | tpep_dropoff_datetime | passenger_count | trip_distance | pickup_longitude | pickup_latitude | |
|---|---|---|---|---|---|---|
| 0 | 2016-05-01 | 2016-05-01 00:17:31 | 1 | 3.60 | -73.985901 | 40.768040 |
| 1 | 2016-05-01 | 2016-05-01 00:07:31 | 1 | 1.68 | -73.991577 | 40.744751 |
| 2 | 2016-05-01 | 2016-05-01 00:07:01 | 6 | 1.09 | -73.993073 | 40.741573 |
| 3 | 2016-05-01 | 2016-05-01 00:19:47 | 1 | 4.21 | -73.991943 | 40.684601 |
| 4 | 2016-05-01 | 2016-05-01 00:06:39 | 1 | 0.56 | -74.005280 | 40.740192 |
data.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 11836853 entries, 0 to 11836852 Data columns (total 19 columns): # Column Dtype --- ------ ----- 0 VendorID int64 1 tpep_pickup_datetime datetime64[ns] 2 tpep_dropoff_datetime datetime64[ns] 3 passenger_count int64 4 trip_distance float64 5 pickup_longitude float64 6 pickup_latitude float64 7 RatecodeID int64 8 store_and_fwd_flag object 9 dropoff_longitude float64 10 dropoff_latitude float64 11 payment_type int64 12 fare_amount float64 13 extra float64 14 mta_tax float64 15 tip_amount float64 16 tolls_amount float64 17 improvement_surcharge float64 18 total_amount float64 dtypes: datetime64[ns](2), float64(12), int64(4), object(1) memory usage: 1.7+ GB
Как видно, файл очень большой, но нас интересуют только данные о месте начала поездки и о времени, поэтому можно считывать не все данные.
В каждом из файлов содержится следующая информация о поездках:
1. время начала поездки
2. время окончания поездки
5,6. долгота и широта точки начала поездки
9,10. долгота и широта точки окончания поездки
3. количество пассажиров
4. расстояние по счётчику
7. тип тарифа (одна из шести категорий)
11. способ оплаты (одна из шести категорий)
9. стоимость поездки по счётчику
12. налог на счётчик
13. доплата за поездки в пиковые часы и ночью
14. доплата за проезд по платным дорогам
17. доплата за проезд, взимаемая с каждой поездки с января 2015
15. размер чаевых
16. общая стоимость поездки
18. Общая сумма, взимаемая с пассажиров. Не включает в себя наличные чаевые.
9. бинарный флаг, показывающий, были ли данные о поездке получены немедленно после её окончания, или какое-то время хранились в памяти автомобиля.
data.describe()
| VendorID | passenger_count | trip_distance | pickup_longitude | pickup_latitude | RatecodeID | dropoff_longitude | dropoff_latitude | payment_type | fare_amount | extra | mta_tax | tip_amount | tolls_amount | improvement_surcharge | total_amount | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| count | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 | 1.183685e+07 |
| mean | 1.529075e+00 | 1.662265e+00 | 6.110069e+00 | -7.304978e+01 | 4.024234e+01 | 1.045014e+00 | -7.311869e+01 | 4.028127e+01 | 1.349840e+00 | 1.334066e+01 | 3.291893e-01 | 4.972577e-01 | 1.831982e+00 | 3.422044e-01 | 2.996838e-01 | 1.664095e+01 |
| std | 4.991539e-01 | 1.304418e+00 | 5.155162e+03 | 8.212003e+00 | 4.523904e+00 | 6.203209e-01 | 7.902544e+00 | 4.353644e+00 | 4.942504e-01 | 1.188799e+01 | 4.585669e-01 | 4.473578e-02 | 2.665043e+00 | 1.760303e+00 | 1.300760e-02 | 1.450379e+01 |
| min | 1.000000e+00 | 0.000000e+00 | 0.000000e+00 | -1.650819e+02 | -4.456935e+01 | 1.000000e+00 | -1.295732e+02 | 0.000000e+00 | 1.000000e+00 | -4.100000e+02 | -3.900000e+01 | -2.700000e+00 | -8.000000e+01 | -1.500000e+01 | -3.000000e-01 | -4.103000e+02 |
| 25% | 1.000000e+00 | 1.000000e+00 | 1.000000e+00 | -7.399179e+01 | 4.073611e+01 | 1.000000e+00 | -7.399116e+01 | 4.073457e+01 | 1.000000e+00 | 6.500000e+00 | 0.000000e+00 | 5.000000e-01 | 0.000000e+00 | 0.000000e+00 | 3.000000e-01 | 8.760000e+00 |
| 50% | 2.000000e+00 | 1.000000e+00 | 1.720000e+00 | -7.398136e+01 | 4.075375e+01 | 1.000000e+00 | -7.397922e+01 | 4.075439e+01 | 1.000000e+00 | 9.500000e+00 | 0.000000e+00 | 5.000000e-01 | 1.260000e+00 | 0.000000e+00 | 3.000000e-01 | 1.230000e+01 |
| 75% | 2.000000e+00 | 2.000000e+00 | 3.270000e+00 | -7.396588e+01 | 4.076862e+01 | 1.000000e+00 | -7.396164e+01 | 4.077012e+01 | 2.000000e+00 | 1.500000e+01 | 5.000000e-01 | 5.000000e-01 | 2.450000e+00 | 0.000000e+00 | 3.000000e-01 | 1.836000e+01 |
| max | 2.000000e+00 | 9.000000e+00 | 1.200000e+07 | 1.184089e+02 | 5.366724e+01 | 9.900000e+01 | 0.000000e+00 | 7.012994e+01 | 5.000000e+00 | 8.007000e+03 | 7.900000e+01 | 4.030000e+01 | 7.500000e+02 | 1.120680e+03 | 3.100000e+00 | 8.007800e+03 |
data.shape
(11836853, 6)
За один месяц имеется почти 12 миллионов записей о поездках.
В основном данные типа float, т.к тут много информации о ценах, расстоянии и координатах.
Считаем тоже самое, но теперь только нужные столбцы.
data = pd.read_csv("yellow_taxi_data\yellow_tripdata_2016-05.csv" , usecols = range(1,7) ,sep=',', header=0,
parse_dates=['tpep_pickup_datetime', 'tpep_dropoff_datetime'], infer_datetime_format=True,
dayfirst=True, engine='c', encoding='cp1251') #считываем только нужные столбцы
data.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 11836853 entries, 0 to 11836852 Data columns (total 6 columns): # Column Dtype --- ------ ----- 0 tpep_pickup_datetime datetime64[ns] 1 tpep_dropoff_datetime datetime64[ns] 2 passenger_count int64 3 trip_distance float64 4 pickup_longitude float64 5 pickup_latitude float64 dtypes: datetime64[ns](2), float64(3), int64(1) memory usage: 541.8 MB
Тем самым, мы сократили количество использованной памяти почти в 4 раза.
Теперь необходимо сделать агрегацию даннных по областям Нью-Йорка.
В данный момент данные имеют вид:
data.head()
| tpep_pickup_datetime | tpep_dropoff_datetime | passenger_count | trip_distance | pickup_longitude | pickup_latitude | |
|---|---|---|---|---|---|---|
| 0 | 2016-05-01 | 2016-05-01 00:17:31 | 1 | 3.60 | -73.985901 | 40.768040 |
| 1 | 2016-05-01 | 2016-05-01 00:07:31 | 1 | 1.68 | -73.991577 | 40.744751 |
| 2 | 2016-05-01 | 2016-05-01 00:07:01 | 6 | 1.09 | -73.993073 | 40.741573 |
| 3 | 2016-05-01 | 2016-05-01 00:19:47 | 1 | 4.21 | -73.991943 | 40.684601 |
| 4 | 2016-05-01 | 2016-05-01 00:06:39 | 1 | 0.56 | -74.005280 | 40.740192 |
New_York_rect = ((-74.25559,-73.70001), (40.49612,40.91553)) #границы Нью-Йорка
Удалим записи которые имеют некоректные данные
data = data.loc[(data['passenger_count']!= 0) & (data['trip_distance']!= 0)] #число пассажиров не 0, длина поездки не 0
data = data.loc[(data['tpep_pickup_datetime']!=data['tpep_dropoff_datetime'])]#время начала не равно времени конца поездки
#поездка в границах Нью-Йорка
data = data.loc[(((data['pickup_longitude'] > New_York_rect[0][0]) & (data['pickup_longitude'] < New_York_rect[0][1]))
&((data['pickup_latitude'] > New_York_rect[1][0]) & (data['pickup_latitude'] < New_York_rect[1][1]) ))]
data.shape # удалили почти 200к записей
(11626521, 6)
data['hour'] = data['tpep_pickup_datetime'].dt.hour #сделаем признаки для дней и времени
data['date'] = data['tpep_pickup_datetime'].dt.date
def cord_to_index(data1 , data2): #функция принимает 2 серии с координатами и возвращает индекс региона
w = (New_York_rect[0][1] - New_York_rect[0][0]) / 50 # ширина области
h = (New_York_rect[1][1] - New_York_rect[1][0]) / 50 # высота области
ind_w = (data1 - New_York_rect[0][0]) // w #номер по ширие т.е ячейка от 0 до 49
ind_h = (data2 - New_York_rect[1][0]) // h #номер по высоте т.е ячейка от 0 до 49
#нумерация в файле regions происходит слево-направо , снизу - вверх ,
return 50*ind_w + ind_h + 1
data['region_id']= cord_to_index(data["pickup_longitude"] , data["pickup_latitude"]).astype('int')
data
| tpep_pickup_datetime | tpep_dropoff_datetime | passenger_count | trip_distance | pickup_longitude | pickup_latitude | hour | date | region_id | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 2016-05-01 00:00:00 | 2016-05-01 00:17:31 | 1 | 3.60 | -73.985901 | 40.768040 | 0 | 2016-05-01 | 1233 |
| 1 | 2016-05-01 00:00:00 | 2016-05-01 00:07:31 | 1 | 1.68 | -73.991577 | 40.744751 | 0 | 2016-05-01 | 1180 |
| 2 | 2016-05-01 00:00:00 | 2016-05-01 00:07:01 | 6 | 1.09 | -73.993073 | 40.741573 | 0 | 2016-05-01 | 1180 |
| 3 | 2016-05-01 00:00:00 | 2016-05-01 00:19:47 | 1 | 4.21 | -73.991943 | 40.684601 | 0 | 2016-05-01 | 1173 |
| 4 | 2016-05-01 00:00:00 | 2016-05-01 00:06:39 | 1 | 0.56 | -74.005280 | 40.740192 | 0 | 2016-05-01 | 1130 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 11836837 | 2016-05-31 23:53:01 | 2016-06-01 00:00:24 | 2 | 1.89 | -73.947876 | 40.774818 | 23 | 2016-05-31 | 1384 |
| 11836838 | 2016-05-31 23:53:02 | 2016-06-01 00:19:36 | 1 | 15.21 | -73.789467 | 40.646942 | 23 | 2016-05-31 | 2068 |
| 11836839 | 2016-05-31 23:53:02 | 2016-05-31 23:59:07 | 1 | 1.80 | -74.008659 | 40.719723 | 23 | 2016-05-31 | 1127 |
| 11836840 | 2016-05-31 23:53:03 | 2016-06-01 00:05:47 | 1 | 2.28 | -73.981644 | 40.758862 | 23 | 2016-05-31 | 1232 |
| 11836841 | 2016-05-19 03:29:16 | 2016-05-19 03:31:24 | 1 | 0.60 | -73.948959 | 40.586658 | 3 | 2016-05-19 | 1361 |
11626521 rows × 9 columns
data['region_id'].max() #есть области где нет поездок
2487
data_count = data.pivot_table(values =['trip_distance'] , index = ['date','hour'] , columns = ['region_id'] , aggfunc = 'count' , fill_value = 0 , margins_name = 'count')
data_count = data_count.rename(columns = {'trip_distance':'count'})
data_count.head(10)
| count | ||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| region_id | 15 | 16 | 25 | 27 | 32 | 45 | 49 | 61 | 77 | 80 | ... | 2476 | 2478 | 2479 | 2481 | 2482 | 2483 | 2484 | 2485 | 2486 | 2487 | |
| date | hour | |||||||||||||||||||||
| 2016-05-01 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
| 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
| 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
| 4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
| 5 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
| 6 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
| 7 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
| 8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
| 9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
10 rows × 1217 columns
Сделаем нормальные индексы
data_count['new_date'] = data_count.index.map(lambda x: pd.datetime.datetime(x[0].year, x[0].month, x[0].day , x[1]))
data_count = data_count.set_index('new_date')
Добавим нулевые столбцы, в которых в этом месяце не было поездок.
# Список регионов с нулевым числом поездок
zeros_list = [('count',y) for y in list(set(range(1,2501)) - set([ x[1] for x in data_count.columns]))]
# Массив нулевых элементов по регоинам и меткам даты и времени
zeros_data = np.zeros(shape=(len(zeros_list), data_count.shape[0]), dtype=int)
# Датафрейм нулевых элементов по регоинам и меткам даты и времени
zeros_dataframe = pd.DataFrame(dict(zip(zeros_list, zeros_data)))
# Переназначение индексов
zeros_dataframe.index=data_count.index
# Расширенный датафрейм количества поездок по всем регионам и меткам даты и времени
data_count = pd.concat([data_count, zeros_dataframe], axis=1, join='inner')
# Сортировка столбцов
data_count = data_count.reindex(sorted(data_count.columns), axis=1)
data_count
| count | |||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| region_id | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | ... | 2491 | 2492 | 2493 | 2494 | 2495 | 2496 | 2497 | 2498 | 2499 | 2500 |
| new_date | |||||||||||||||||||||
| 2016-05-01 00:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-01 01:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-01 02:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-01 03:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-01 04:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2016-05-31 19:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 20:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 21:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 22:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 23:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
744 rows × 2500 columns
Теперь можно сохранить данные за месяц.
data_count.to_csv('timed_data//timed_yellow_tripdata_2016-05.csv' , index = True , header = False)
data_count
| count | |||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| region_id | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | ... | 2491 | 2492 | 2493 | 2494 | 2495 | 2496 | 2497 | 2498 | 2499 | 2500 |
| new_date | |||||||||||||||||||||
| 2016-05-01 00:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-01 01:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-01 02:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-01 03:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-01 04:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2016-05-31 19:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 20:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 21:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 22:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 23:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
744 rows × 2500 columns
Чтобы было удобно обрабатывать файлы разных месяцев, оьбъединим все в 1 общую функцию.
def prepare_data(filename , src_dir, des_dir):
data = pd.read_csv(os.path.join(src_dir,filename) , usecols = range(1,7),sep=',', header=0,
parse_dates=['tpep_pickup_datetime', 'tpep_dropoff_datetime'], infer_datetime_format=True,
dayfirst=True, engine='c', encoding='cp1251') #считываем только нужные столбцы
New_York_rect = ((-74.25559,-73.70001), (40.49612,40.91553))
#число пассажиров не 0, длина поездки не 0
data = data.loc[(data['passenger_count']!= 0) & (data['trip_distance']!= 0)]
#время начала не равно времени конца поездки
data = data.loc[(data['tpep_pickup_datetime']!=data['tpep_dropoff_datetime'])]
#поездка в границах Нью-Йорка
data = data.loc[(((data['pickup_longitude'] > New_York_rect[0][0]) &
(data['pickup_longitude'] < New_York_rect[0][1]))&
((data['pickup_latitude'] > New_York_rect[1][0]) &
(data['pickup_latitude'] < New_York_rect[1][1]) ))]
data['hour'] = data['tpep_pickup_datetime'].dt.hour
data['date'] = data['tpep_pickup_datetime'].dt.date
data['region_id']= cord_to_index(data["pickup_longitude"] , data["pickup_latitude"]).astype('int')
data_count = data.pivot_table(values =['trip_distance'] , index = ['date','hour'] , columns = ['region_id'] , \
aggfunc = 'count' , fill_value = 0)
data_count = data_count.rename(columns = {'trip_distance':'count'})
data_count['new_date'] = data_count.index.map(lambda x: pd.datetime.datetime(x[0].year, x[0].month, x[0].day , x[1]))
data_count = data_count.set_index('new_date')
# Список регионов с нулевым числом поездок
zeros_list = [('count',y) for y in list(set(range(1,2501)) - set([ x[1] for x in data_count.columns]))]
# Массив нулевых элементов по регоинам и меткам даты и времени
zeros_data = np.zeros(shape=(len(zeros_list), data_count.shape[0]), dtype=int)
# Датафрейм нулевых элементов по регоинам и меткам даты и времени
zeros_dataframe = pd.DataFrame(dict(zip(zeros_list, zeros_data)))
# Переназначение индексов
zeros_dataframe.index=data_count.index
# Расширенный датафрейм количества поездок по всем регионам и меткам даты и времени
data_count = pd.concat([data_count, zeros_dataframe], axis=1, join='inner')
data_count = data_count.reindex(sorted(data_count.columns), axis=1)
data_count.to_csv(os.path.join(des_dir ,'timed_' + filename) , index = True , header = False)
Обработаем еще несколько месяцев
import tqdm
src_dir = 'yellow_taxi_data'
des_dir = 'timed_data'
for filename in tqdm.tqdm(os.listdir(src_dir)[2:]):
prepare_data(filename , src_dir, des_dir)
100%|████████████████████████████████████████████████████████████████████████████████████| 6/6 [06:15<00:00, 62.54s/it]
Подключим нужные библиотеки
import os
import matplotlib.pyplot as plt
import matplotlib
import conda
import folium
import seaborn as sns
import matplotlib.ticker as ticker
import matplotlib , branca
from matplotlib import cm
from matplotlib.colors import ListedColormap
from folium import Choropleth, Circle, Marker
from folium.plugins import HeatMap, MarkerCluster
plt.style.use('ggplot')
%matplotlib inline
conda_file_dir = conda.__file__
conda_dir = conda_file_dir.split('lib')[0]
proj_lib = os.path.join(os.path.join(conda_dir, 'share'), 'proj')
os.environ["PROJ_LIB"] = proj_lib
from mpl_toolkits.basemap import Basemap
import warnings #отключаем предупреждения
warnings.filterwarnings("ignore")
Будем использовать 4 месяца. С 2016-02 по 2016-05.
datas = []
for filename in os.listdir(des_dir)[1:5]:
data = pd.read_csv(os.path.join(des_dir,filename), header = None,
index_col = 0,parse_dates=[0], infer_datetime_format=True,
engine='c', encoding='cp1251')
datas.append(data)
full_data = pd.concat(datas)
full_data
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | ... | 2491 | 2492 | 2493 | 2494 | 2495 | 2496 | 2497 | 2498 | 2499 | 2500 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | |||||||||||||||||||||
| 2016-02-01 00:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-02-01 01:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-02-01 02:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-02-01 03:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-02-01 04:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2016-05-31 19:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 20:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 21:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 22:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 2016-05-31 23:00:00 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2904 rows × 2500 columns
Для начала отрисуем график числа поездок из выбранной нами области.
fig, ax = plt.subplots(figsize=(20,8))
ax.plot(range(full_data.shape[0]) , full_data[1231])
ax.set_xlabel(xlabel='День и час посадки в такси')
ax.set_ylabel(ylabel='Число поездок')
ax.set_title('График количества поездок жёлтого такси за февраль-май 2016 в одной из центральных ячеек Нью-Йорка')
ax.xaxis.set_major_locator(ticker.MultipleLocator(168))
ax.xaxis.set_minor_locator(ticker.MultipleLocator(24))
ax.yaxis.set_major_locator(ticker.MultipleLocator(100))
ax.yaxis.set_minor_locator(ticker.MultipleLocator(50))
ax.grid(which='major', axis='both', color = 'gray', linewidth='0.5', linestyle='-')
ax.grid(which='minor', axis='both', color = 'gray', linewidth='0.5', linestyle=':')
ax.set_xlim(0, full_data.shape[0])
ax.set_xticklabels([0] + [x.day for x in full_data.index.date][::168], minor=False, fontdict={'fontsize':12})
ax.tick_params(axis='x', which='minor', pad=0)
ax.tick_params(axis='x', which='major', pad=15)
fig.autofmt_xdate(bottom=0.2, rotation=0, which='minor')
fig.autofmt_xdate(bottom=0.2, rotation=0, which='major')
plt.show()
Вывод : из этого графика явно видно выраженную дневную и недельную сезонность.
Это логично, люди чаще пользуются такси в вечерное время (как и в будние дни).
Поскольку ряд имеет сезонность, он не является стационарным.
Визуально в течении первых 3х месяцев среднее значение и дисперсия ряда выглядят постоянно.
Но вот в мае, значения становятся чуть меньшье (также как и дисперсия).
Возможно это связано с тем, что в мае уже достаточно тепло, чтобы ходить пешком).
Скорее всего это следствие годовой сезонности, которую мы не можем рассмотреть на отрезке в несколько месяцев.
Появляется отличная возможнось с помощью статистических критериев проверить гипотезу о равенстве среднего числа поездок за май месяц и числа поездок за предыдущие месяцы.
Посмотрим какое распределение имеет число поездок за весь этот период.
full_data[1231].hist(bins=40,figsize=(20,8), alpha=0.7)
plt.axvline(x=full_data[1231].mean(), c='g', label='Mean ' + str(round(full_data[1231].mean())), linestyle='dashed')
plt.axvline(x=full_data[1231].median(), c='y', label='Median' + str(round(full_data[1231].median())), linestyle='dashed')
plt.legend()
<matplotlib.legend.Legend at 0x1b06ee85748>
Видно, что оно имеет несколько пиков, а точнее 2 основных пика. Основная масса приходится на большие значения (около 1000).
Распределение совсем не похоже на нормальное.
Среднее сильно отличается от мод.
Можем сказать, что чаще всего люди используют такси либо "мало", либо "много".
Теперь взгляним на картину в общем и поиграемся с картами.
Отрисуем суммарное количество поездок за май.
sum_data = np.sum(full_data , axis = 0)#суммируем число поездок
sum_data= (np.array(sum_data).reshape((50,50))).T
#широты и высоты областей
lx = np.linspace(New_York_rect[0][0] ,New_York_rect[0][1],51)
ly = np.linspace(New_York_rect[1][0] ,New_York_rect[1][1],51)
#сетка в координатах Земли
xx , yy = np.meshgrid(lx,ly)
center = [np.mean(cor) for cor in New_York_rect] #центр
fig = plt.figure(num=None, figsize=(15, 15))
#создаем карту
map_basemap = Basemap(llcrnrlat=(New_York_rect[1][0]),
llcrnrlon=(New_York_rect[0][0]),
urcrnrlat=(New_York_rect[1][1]),
urcrnrlon=(New_York_rect[0][1]),
resolution='f',
projection='lcc',
area_thresh=5000,
lat_0=center[0],
lon_0=center[1],
epsg=3857
)
my_cmap = plt.get_cmap('viridis' , 256)
newcolors = my_cmap(np.linspace(0, 1, 256))
pink = np.array([248/256, 24/256, 148/256, 1])
newcolors[:1, :] = pink
newcmp = ListedColormap(newcolors)
#сетка в координатах графика
llx, lly = map_basemap(xx, yy)
colormesh = map_basemap.pcolormesh(llx,lly, data = sum_data,cmap = newcmp , alpha = 0.3 )
cb = map_basemap.colorbar(colormesh, location='bottom', label="Summary taxi ride")
#выбираем текстуру города
map_basemap.arcgisimage(service='World_Street_Map', xpixels=1500)
#отрисуем границы областей
map_basemap.drawparallels(np.linspace(New_York_rect[1][0] ,New_York_rect[1][1],51) ,linewidth=0.25, color='0.5')
map_basemap.drawmeridians(np.linspace(New_York_rect[0][0] ,New_York_rect[0][1],51) , linewidth=0.25, color='0.5')
plt.title("Static NY Map")
plt.show()
Отрисуем среднее в час количество поездок за май.
mean_data = np.array(np.mean(full_data , axis = 0)).reshape((50,50)).T #находим среднее число поездок
def get_geojson_grid(upper_right, lower_left, mean_matrix, n=50):
#Создает сетку в формате GeoJson, чтобы дальше отрисовать ее с помощью folium
#Аргументы : Правый верхний угол ,Нижний левый, Матрица со средним числом поездок,на сколько частей разбивать
geo_json = {"type": "FeatureCollection",
"features":[]}
lat_steps = np.linspace(lower_left[0], upper_right[0], n+1)
lon_steps = np.linspace(lower_left[1], upper_right[1], n+1)
lat_stride = lat_steps[1] - lat_steps[0]
lon_stride = lon_steps[1] - lon_steps[0]
for i , lat in enumerate(lat_steps[:-1]):
for k , lon in enumerate(lon_steps[:-1]):
# Define dimensions of box in grid
upper_left = [lon, lat + lat_stride]
upper_right = [lon + lon_stride, lat + lat_stride]
lower_right = [lon + lon_stride, lat]
lower_left = [lon, lat]
# Define json coordinates for polygon
coordinates = [
upper_left,
upper_right,
lower_right,
lower_left,
upper_left
]
grid_feature = {
"type":"Feature",
"geometry":{
"type":"Polygon",
"coordinates": [coordinates],
},
"properties":{
"mean": mean_matrix[i][k]
},
}
geo_json["features"].append(grid_feature)
return geo_json
#создаем палитру для отрисовки
colormap = branca.colormap.LinearColormap(colors=['white','yellow' ,'green','red'], index=[0,5,50,np.max(mean_data)],
vmin = 0 ,vmax = np.max(mean_data))
colormap
Начиная с 5 поездок среднее будет показываться желтым цветом , c 50 - зеленым.
upper_right = [New_York_rect[1][1] , New_York_rect[0][1]]
lower_left = [New_York_rect[1][0] , New_York_rect[0][0]]
grid = get_geojson_grid(upper_right,lower_left, mean_data, n=50)
gj = folium.GeoJson(grid, style_function=lambda x: {#'fillColor': colorscale(x["properties"]["mean"]),
'color':colormap(x["properties"]["mean"]),
'weight': 1,
'dashArray': '5, 5',
'opacity': 0.6,
'fillOpacity': 0.5,
}, name="Mean Taxi Ride")
# Create a map
m_NY = folium.Map(location=center[::-1], tiles='openstreetmap', zoom_start=10)
our_cell = [40.748391, -73.985543] #координаты нашей ячейки
Marker(our_cell).add_to(m_NY) #Добавляем маркер нашей ячейки
<folium.map.Marker at 0x1e279d571d0>
m_NY.add_child(gj)
Получили красивую карту.
Судя по карте основная масса поездок приходится на Манхэттен. Что логично, т.к это деловой центр города.
Но также достаточное количество поездок приходится на аэропорты.
Ячеек из неадекватных мест(водоемов и т.д) вроде бы нет, так что скорее всего мы все сделали правильно.
Теперь разложим ряд на компоненты и отрисуем STL декомпозицию.
import statsmodels.api as sm
from pylab import rcParams
rcParams['figure.figsize'] = 12, 7
sm.tsa.seasonal_decompose(full_data[1231]).plot()
plt.show()
Вывод:
Визуально ряд не имеет ярко выраженного тренда. Значения колеблются в одних и тех же пределах.
Но в майский месяц значения более сжаты, будто дисперсия рядо уменьшилась.
Сезонная компонента состоит из дневной сезонности.
Недельняя сезонность ушла в тренд. Видно как тренд раскладывается на отрезки, где минимумы приходятся на выходные дни.
from scipy import stats
Возьмем функцию с лекции.
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries, window = 24, cutoff = 0.05):
#Determing rolling statistics
rolmean = timeseries.rolling(window).mean()
rolstd = timeseries.rolling(window).std()
#Plot rolling statistics:
fig = plt.figure(figsize=(12, 4))
orig = plt.plot(timeseries, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show()
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.show()
#Perform Dickey-Fuller test:
print('Results of Dickey-Fuller Test:')
dftest = adfuller(timeseries.values,autolag='AIC' )
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
pvalue = dftest[1]
if pvalue < cutoff:
print('p-value = %.4f. The series is likely stationary.' % pvalue)
else:
print('p-value = %.4f. The series is likely non-stationary.' % pvalue)
print(dfoutput)
test_stationarity(full_data[1231])
Results of Dickey-Fuller Test: p-value = 0.0000. The series is likely stationary. Test Statistic -6.794182e+00 p-value 2.323501e-09 #Lags Used 2.800000e+01 Number of Observations Used 2.875000e+03 Critical Value (1%) -3.432627e+00 Critical Value (5%) -2.862546e+00 Critical Value (10%) -2.567305e+00 dtype: float64
Исходя из критерия Дики-Фулера, получается что ряд стационарен (т.к значение p-value 2.323501e-09).
Но мы знаем, что это не так, т.к есть сезонность.
Возможно, он не сработал из-за большой зашумленности данных.
Визуально, очевидно, что данные распределены не нормально, но давайте в этом убедимся с помощью критерия Шапиро-Уилка:
$H_0\colon$ данные распределны нормально.
$H_1\colon$ не нормально.
stats.shapiro(full_data[1231])
ShapiroResult(statistic=0.9452146887779236, pvalue=1.4312167991741817e-31)
Критерий уверенно отвергает нулевую гипотезу в пользу альтернативу, с достигаемым уровнем значимости pvalue=1.4312167991741817e-31.
Мы убедились, что данные распределны не нормально.
Сначала просто посмотрим на средние.
for i in range(2,6):
print(f'Mean Count of Ride in {i} month: {full_data[1231][full_data.index.month==i].mean()}')
print('Full mean:',full_data[1231].mean())
Mean Count of Ride in 2 month: 722.125 Mean Count of Ride in 3 month: 740.8494623655914 Mean Count of Ride in 4 month: 700.8013888888889 Mean Count of Ride in 5 month: 657.9153225806451 Full mean: 705.1849173553719
Видно что в мае среднее прилично меньше.
full_data[1231][full_data.index.month==5].std(), full_data[1231][full_data.index.month!=5].std()
(320.7736435098047, 364.34509164835004)
Дисперсии также отличаются.
Далее, будем проверять гипотезу о равенсте среднего числа поездок против двусторонней альтернативы. (Вместо того, чтобы брать одностороннюю альтернативу).
Это делается для того, чтобы мы не переобучалсь, т.к мы уже видели данные и знаем, что среднее скорее всего меньше.
Попробуем использовать Критерий Стьюдента, для проверки гипотезы о том, что среднее количество поездок за май не равно среднему количеству поездок за предыдущие месяцы.
H0: число заказов в мае равно числу заказов за 2-4 месяц.
H1: число заказов в мае не равно числу заказов за 2-4 месяц.
Хотя он и требует нормальности данных и равенства дисперсий, посмотрим на результаты.
t_stat = stats.ttest_ind(full_data[1231][full_data.index.month==5], full_data[1231][full_data.index.month!=5])
t_stat
Ttest_indResult(statistic=-4.226719605985014, pvalue=2.444402377784712e-05)
Критерий уверенно отвергает нулевую гипотезу в пользу альтернативу, с достигаемым уровнем значимости pvalue=2.444402377784712e-05.
Критерий говорит, что средние значения все таки различны.
Но чтобы быть уверенными в результатах, используем критерий, который не требует никаких дполнительных условий от данных.
Воспользуем одним из непараметрических критериев - Критерием Манна-Уитни.
stats.mannwhitneyu(full_data[1231][full_data.index.month==5], full_data[1231][full_data.index.month!=5])
MannwhitneyuResult(statistic=701769.5, pvalue=1.2434458841687715e-07)
Критерий также уверенно отвергает нулевую гипотезу, с достигаемым уровнем значимости pvalue=1.2434458841687715e-07.
Можем сделать вывод: число поездок в мае значительно уменьшилось.
from statsmodels.sandbox.stats.multicomp import multipletests
Проведем попарные тесты на равенство средних.
stat_data = []
for i in range(2,6):
for j in range(2,6):
if i >= j:
continue
stat, p = stats.mannwhitneyu(full_data[1231][full_data.index.month==i], full_data[1231][full_data.index.month==j])
stat_data.append([i, j, p])
stat_data = pd.DataFrame(stat_data, columns=['month_1', 'month_2', 'p_val'])
stat_data
| month_1 | month_2 | p_val | |
|---|---|---|---|
| 0 | 2 | 3 | 1.942510e-01 |
| 1 | 2 | 4 | 4.317316e-02 |
| 2 | 2 | 5 | 7.558823e-06 |
| 3 | 3 | 4 | 4.385482e-03 |
| 4 | 3 | 5 | 3.848086e-08 |
| 5 | 4 | 5 | 2.150166e-03 |
С помощью Критерия Манна-Уитни посчитали p-value, для гипотез о равенстве средних.
Теперь применим Метод Холма для поправки на множественную проверку гипотез.
reject, p_corrected, a1, a2 = multipletests(stat_data.p_val,
alpha = 0.05,
method = 'holm')
stat_data['p_corrected'] = p_corrected
stat_data['reject'] = reject
stat_data
| month_1 | month_2 | p_val | p_corrected | reject | |
|---|---|---|---|---|---|
| 0 | 2 | 3 | 1.942510e-01 | 1.942510e-01 | False |
| 1 | 2 | 4 | 4.317316e-02 | 8.634632e-02 | False |
| 2 | 2 | 5 | 7.558823e-06 | 3.779411e-05 | True |
| 3 | 3 | 4 | 4.385482e-03 | 1.315645e-02 | True |
| 4 | 3 | 5 | 3.848086e-08 | 2.308852e-07 | True |
| 5 | 4 | 5 | 2.150166e-03 | 8.600663e-03 | True |
Как видно, гипотезы о равенстве средних отклоняются, для случаев кроме 2-3, 2-4 месяца.
В остальных же случаях считаем, что средние не равны.
Интересно, что средние 2-3, 2-4 равны, а 3-4 нет.
Будем делать предсказания спроса на такси на ближайшие сутки.
from sklearn.metrics import mean_squared_error
from math import sqrt
def rmse(y_actual, y_predicted):
rmse = sqrt(mean_squared_error(y_actual, y_predicted))
return rmse
# Относительная ошибка - хорошая метрика для бизнеса
def mean_absolute_percentage_error(y_true, y_pred):
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
from statsmodels.tsa.api import ExponentialSmoothing, SimpleExpSmoothing, Holt
from sklearn.model_selection import TimeSeriesSplit
Наивный baseline
w_hours = 24
train = full_data[1231].iloc[:-w_hours]
val = full_data[1231].iloc[-w_hours:]
preds = [train.mean()] * w_hours #просто среднее
print('RMSE - ', rmse(val, preds))
print('MAPE - ', mean_absolute_percentage_error(val, preds))
RMSE - 366.0676621035622 MAPE - 199.22765781550595
preds = train.iloc[-24:].values #то же, что и вчера
print('RMSE - ', rmse(val, preds))
print('MAPE - ', mean_absolute_percentage_error(val, preds))
RMSE - 297.4780804249841 MAPE - 61.073151812695606
Видим довольно большие ошибки. Попробуем что-нибудь посильнее.
Посчитаем качество тройного экспоненциального сглаживания на кросс-валидации и примем его за baseline.
fit1 = ExponentialSmoothing(np.asarray(train) , seasonal_periods=24, seasonal='add').fit()
forecast = pd.Series(fit1.forecast(len(val)))
forecast.index = val.index
plt.plot(val)
plt.plot(forecast)
[<matplotlib.lines.Line2D at 0x1e26a2a09b0>]
print('RMSE - ', rmse(val, forecast))
print('MAPE - ', mean_absolute_percentage_error(val, forecast))
RMSE - 479.3803030806802 MAPE - 128.56354218180576
Предсказания вообще выглядят довольно неадекватно, принимая значения даже меньше 0. Но посмотрим какое качество будет в среднем.
errors_rmse = []
errors_mape = []
tscv = TimeSeriesSplit(n_splits=10)
for train_idx, test_idx in tscv.split(full_data[1231]):
#print('len_train - ', len(train_idx))
#print('len_test - ', len(test_idx))
model = ExponentialSmoothing(full_data[1231][train_idx] ,seasonal_periods=24 , seasonal='add').fit()
forecast = pd.Series(model.forecast(len(test_idx[:w_hours])))
actual = full_data[1231][test_idx[:w_hours]] #из тестовой выборки берем всего лишь сутки,
#т.к в итоге хотим предсказывать имеено такию длительность.
errors_rmse.append(rmse(actual.values, forecast.values))
errors_mape.append(mean_absolute_percentage_error(actual.values, forecast.values))
np.mean(errors_mape), np.mean(errors_rmse)
(50.951444229587125, 247.88907037012495)
rmse меньше чем было раньше, посмотрим что можно с этим делать.
from sklearn.linear_model import LinearRegression, Ridge
Идея: давайте построим обычную линейную модель на придуманных нами регрессионных признаках. А информацию, которая останется в остатках от регресси, будем предсказывать с помощью ARIM'ы. Конечно, можно передать регрессионные признаки в Arima (c помощью параметра exog), но Arima и так учится слишком большое количество времени, а прирост качества должен из-за этого быть очень незначительный.
В рамках моделей ARIMA можно учесть только одну из сезонностей.
Обычно в таких случаях сезонность с самым маленьким периодом явно моделируют с помощью аримы, а все остальные учитывают за счёт регрессионной компоненты.
Недельную сезонность будем учитывать с помощью гармонических признаков.
from category_encoders.target_encoder import TargetEncoder
from sklearn.preprocessing import StandardScaler
def regression_matrix(data , T , K, lag_count, train_ind, test_ind): #матрица регрессионных признаков
y_copy = pd.DataFrame(data)
#делим выборку на 2 части
y_copy.loc[y_copy.index[train_ind],'label'] = 'train'
y_copy.loc[y_copy.index[test_ind],'label'] = 'test'
matr_sin = np.zeros((T , K))
matr_cos = np.zeros((T , K))
for t in range(T):#гармонические признаки
for i in range(K):
matr_sin[t,i] = np.sin((t)*2*np.pi*i/168)
matr_cos[t,i] = np.cos((t)*2*np.pi*i/168)
#добавим индикатор дня недели
day_of_week = np.zeros((T , 1))
for day_num in range(len(y_copy.index)): #нам важен только порядок дней недели
day_of_week[k] = day_num % 7
y_copy['day_of_week'] = day_of_week
#добавим индикатор часа.
hour_id_day = np.zeros((T , 1))
for k , index in enumerate(y_copy.index):
hour_id_day[k] = index.hour
y_copy['hour_id_day'] = hour_id_day
# лаговые фичи
lag_period = 24
features = []
for period_mult in range(1, lag_count, 1):
y_copy["lag_period_{}".format(period_mult)] = y_copy[1231].shift(period_mult * lag_period)
features.append("lag_period_{}".format(period_mult))
#лаговые статистики
y_copy['lagf_mean'] = y_copy[features].mean(axis = 1)
target_encoder = TargetEncoder(cols = ['day_of_week','hour_id_day'], #сделаем таргет кодирование для дня недели и часа
return_df = True)
coding_train = target_encoder.fit_transform(y_copy[y_copy.label == 'train'].drop([1231, 'label'], axis = 1),
y_copy[y_copy.label == 'train'][1231])
coding_test = target_encoder.transform(y_copy[y_copy.label == 'test'].drop([1231, 'label'], axis = 1))
scaler = StandardScaler() #сделаем нормализацию данных для линейной модели
coding_train_scaler = scaler.fit_transform(coding_train)
coding_test_scaler = scaler.transform(coding_test)
#объединим все в одну матрицу
return np.hstack((matr_sin, matr_cos, np.vstack((coding_train_scaler, coding_test_scaler))))
Теперь сделаем небольшую кросс-валидацию, для того чтобы подобрать параметры
K - для признаков $ sin([1,…,T] ∗ 2πi/168), \quad cos([1,…,T]∗2πi/168), i \in 1,...,K $, учитывающих недельную сезонность
log_count - число лаговых фичей
def cross_validate_linear(K, lag_count):
errors_rmse = []
errors_mape = []
tscv = TimeSeriesSplit(n_splits=10)
for train_idx, test_idx in tscv.split(full_data[1231]):
len_ = len(train_idx) + len(test_idx)
regression_feat = regression_matrix(full_data[1231][:len_], len_, K, lag_count, train_idx, test_idx)
model = Ridge().fit(regression_feat[train_idx][lag_count * 24:],
full_data[1231][train_idx][lag_count * 24:]) #тут нужно обрезать строки с nan'ми
forecast = model.predict(regression_feat[test_idx[:w_hours]]) #качество считаем не на всей тестовой выборке
actual = full_data[1231][test_idx[:w_hours]] #а на 24 часах, т.к это цель задания
errors_rmse.append(rmse(actual.values, forecast))
errors_mape.append(mean_absolute_percentage_error(actual.values, forecast))
print(f'Params: K - {K}, lag_count - {lag_count}')
print('Score:', np.mean(errors_mape[2:]), np.mean(errors_rmse[2:])) #тут берем с 2, т.к первые обучающие выборки
#совсем маленькие из-за убирания nan'oв
from sklearn.model_selection import ParameterGrid
params = ParameterGrid({'K': [4, 6, 8, 10, 12], 'lag_count': [6, 8, 10]})
for param in params:
cross_validate_linear(*param.values())
Params: K - 4, lag_count - 6 Score: 23.53994949743754 126.5411026775684 Params: K - 4, lag_count - 8 Score: 13.577301735143065 82.13693606046408 Params: K - 4, lag_count - 10 Score: 13.288270959558732 82.38615649601738 Params: K - 6, lag_count - 6 Score: 19.324849845641232 117.36421408362536 Params: K - 6, lag_count - 8 Score: 12.566559810079541 79.05315274683944 Params: K - 6, lag_count - 10 Score: 12.268166731110279 80.92905058816396 Params: K - 8, lag_count - 6 Score: 19.1243358954262 116.43815933908186 Params: K - 8, lag_count - 8 Score: 12.76275593634327 80.75915729018604 Params: K - 8, lag_count - 10 Score: 12.375740796755583 81.91584450083727 Params: K - 10, lag_count - 6 Score: 17.946244635459294 110.43935508492929 Params: K - 10, lag_count - 8 Score: 12.557753790976477 80.82720016025164 Params: K - 10, lag_count - 10 Score: 12.597406533139516 82.37623315835833 Params: K - 12, lag_count - 6 Score: 17.062664385100494 108.45061331860413 Params: K - 12, lag_count - 8 Score: 11.880454673188066 80.4179657335404 Params: K - 12, lag_count - 10 Score: 11.756565787187741 82.54741429477735
plt.plot(actual.values, label = 'y')
plt.plot(forecast, label = 'prediction')
plt.legend()
<matplotlib.legend.Legend at 0x1e27a3b8400>
Теперь график выглядит очень похоже. Параметры возьмем K = 12, lag_count = 8.
Обучим такую модель на 3 месяцах и посмотрим остатки.
K = 12
lag_count = 8
train_idx = full_data.index.month != 5
test_idx = full_data.index.month == 5
regression_feat = regression_matrix(full_data[1231], T, K, lag_count, train_idx, test_idx)
model = Ridge().fit(regression_feat[train_idx][lag_count * 24:],
full_data[1231][train_idx][lag_count * 24:])
y_predict = model.predict(regression_feat[train_idx][lag_count * 24:])
Отрисуем график для трейна.
full_data[1231][train_idx][lag_count * 24:].plot(figsize = (15, 10), label = 'Изначальный ряд')
(y_residual + full_data[1231][train_idx][lag_count * 24:]).plot(label = 'Построенная модель')
plt.legend()
<matplotlib.legend.Legend at 0x1e27a4c42b0>
y_residual = y_predict - full_data[1231][train_idx][lag_count * 24:]
y_residual.plot(figsize = (12, 8))
plt.title('Остатки от линейной модели')
Text(0.5, 1.0, 'Остатки от линейной модели')
На самом деле остатки сильно напоминают белый шум. Явного тренда нет, сезонности тоже.
Проверим несмещенность остатков, с помощью критерия Стьюдента и стационарность с помощью критерия Дики-Фулера.
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(y_residual, 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(y_residual)[1])
Критерий Стьюдента: p=1.000000 Критерий Дики-Фуллера: p=0.000000
Как мы видим, остатки являются несмещенными и стационарными. Критерии уверенно это признают.
Вообще, на этом можно было бы остановиться, но давайте посмотрим, можно ли что-нибудь сделать с помощью Arim'ы
Отрисуем корролеллограммы.
y_residual = pd.DataFrame(y_residual)
rcParams['figure.figsize'] = 12, 7
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(y_residual[1231], lags=48, ax=ax)
plt.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(y_residual[1231], lags=48, ax=ax)
plt.show()
Из графика автокорреляции все таки видно, что некоторая сезонность осталась в данных. Попробуем его продифференцировать.
y_residual['diff_24'] = y_residual[1231] - y_residual[1231].shift(24) #Первая сезонная производная
y_residual['diff_1'] = y_residual[1231] - y_residual[1231].shift(1) #Первая обычная производная
y_residual.head()
| 1231 | diff_24 | diff_1 | |
|---|---|---|---|
| 0 | |||
| 2016-02-09 00:00:00 | 118.315309 | NaN | NaN |
| 2016-02-09 01:00:00 | 49.352546 | NaN | -68.962763 |
| 2016-02-09 02:00:00 | 49.514364 | NaN | 0.161818 |
| 2016-02-09 03:00:00 | 12.892994 | NaN | -36.621370 |
| 2016-02-09 04:00:00 | 30.095321 | NaN | 17.202327 |
Еще раз отрисуем королеллограмму.
Для сезонной производной:
rcParams['figure.figsize'] = 12, 7
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(y_residual['diff_24'][24:], lags=48, ax=ax)
plt.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(y_residual['diff_24'][24:], lags=48, ax=ax)
plt.show()
Для обычной производной:
rcParams['figure.figsize'] = 12, 7
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(y_residual['diff_1'][1:], lags=48, ax=ax)
plt.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(y_residual['diff_1'][1:], lags=48, ax=ax)
plt.show()
Остановимся на обычной проиводной, т.к королеллограмма выглядит более независимой, чем для сезонного случая.
Из первого графика можно сделать вывод что:
Q = 1
q= 1 или 2
Из второго графика:
P = 1
p - тоже возьмем не большим (попробуем 0 - 2), хотя значимых значений намного больше.
И параметры дифференцирования:
d = 1
D = 0
ps = range(0, 3)
d = 1
qs = range(0, 3)
Ps = range(0,2)
D = 0
Qs = range(0,2)
from itertools import product
parameters = product(ps, qs, Ps,Qs)
parameters_list = list(parameters)
len(parameters_list)
36
Подберем параметры Arim'ы
%%time
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in parameters_list:
try:
model=sm.tsa.statespace.SARIMAX(y_residual[1231], order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 24)).fit(disp=-1)
except:
print('wrong parameters:', param)
continue
print('Ready:', param)
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
Ready: (0, 0, 0, 0) Ready: (0, 0, 0, 1) Ready: (0, 0, 1, 0) Ready: (0, 0, 1, 1) Ready: (0, 1, 0, 0) Ready: (0, 1, 0, 1) Ready: (0, 1, 1, 0) Ready: (0, 1, 1, 1) Ready: (0, 2, 0, 0) Ready: (0, 2, 0, 1) Ready: (0, 2, 1, 0) Ready: (0, 2, 1, 1) Ready: (1, 0, 0, 0) Ready: (1, 0, 0, 1) Ready: (1, 0, 1, 0) Ready: (1, 0, 1, 1) Ready: (1, 1, 0, 0) Ready: (1, 1, 0, 1) Ready: (1, 1, 1, 0) Ready: (1, 1, 1, 1) Ready: (1, 2, 0, 0) Ready: (1, 2, 0, 1) Ready: (1, 2, 1, 0) Ready: (1, 2, 1, 1) Ready: (2, 0, 0, 0) Ready: (2, 0, 0, 1) Ready: (2, 0, 1, 0) Ready: (2, 0, 1, 1) Ready: (2, 1, 0, 0) Ready: (2, 1, 0, 1) Ready: (2, 1, 1, 0) Ready: (2, 1, 1, 1) Ready: (2, 2, 0, 0) Ready: (2, 2, 0, 1) Ready: (2, 2, 1, 0) Ready: (2, 2, 1, 1) Wall time: 2min 27s
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head())
parameters aic 18 (1, 1, 1, 0) 22071.578443 19 (1, 1, 1, 1) 22072.582149 30 (2, 1, 1, 0) 22073.566144 22 (1, 2, 1, 0) 22073.566168 23 (1, 2, 1, 1) 22074.560002
print(best_model.summary())
SARIMAX Results
===========================================================================================
Dep. Variable: 1231 No. Observations: 1968
Model: SARIMAX(1, 1, 1)x(1, 0, [], 24) Log Likelihood -11031.789
Date: Sat, 12 Dec 2020 AIC 22071.578
Time: 02:38:16 BIC 22093.916
Sample: 02-09-2016 HQIC 22079.787
- 04-30-2016
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.6447 0.011 56.230 0.000 0.622 0.667
ma.L1 -0.9960 0.002 -405.975 0.000 -1.001 -0.991
ar.S.L24 -0.1356 0.021 -6.517 0.000 -0.176 -0.095
sigma2 4345.8318 90.415 48.065 0.000 4168.621 4523.043
===================================================================================
Ljung-Box (Q): 70.32 Jarque-Bera (JB): 886.21
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 0.81 Skew: -0.10
Prob(H) (two-sided): 0.01 Kurtosis: 6.28
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
В итоге лучшый набор параметров оказался
p = 1, q = 1, P = 1, Q = 0
Отрисуем финальные остатки.
rcParams['figure.figsize'] = 15, 8
plt.subplot(211)
best_model.resid[1:].plot()
plt.ylabel(u'Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(best_model.resid[1:].values.squeeze() , lags=48, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(best_model.resid[1:].values.squeeze() , 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(best_model.resid[1:])[1])
Критерий Стьюдента: p=0.438464 Критерий Дики-Фуллера: p=0.000000
Теперь это точно белый шум. Они также не смещены и стационарны.
Давайте совместим линейную модель и Arim'у и проверим качество на кросс-валидации.
K = 12
lag_count = 8
p = 1
q = 1
P = 1
Q = 0
d = 1
D = 0
errors_rmse = []
errors_mape = []
tscv = TimeSeriesSplit(n_splits=10)
for train_idx, test_idx in tscv.split(full_data[1231]):
len_ = len(train_idx) + len(test_idx)
regression_feat = regression_matrix(full_data[1231][:len_], len_, K, lag_count, train_idx, test_idx)
#обучаем линейную модель
linear_model = Ridge().fit(regression_feat[train_idx][lag_count * 24:],
full_data[1231][train_idx][lag_count * 24:]) #тут нужно обрезать строки с nan'ми
#считаем невязку
residual = linear_model.predict(regression_feat[train_idx][lag_count * 24:]) - \
full_data[1231][train_idx][lag_count * 24:]
#обучам Arim'у
arima_model = sm.tsa.statespace.SARIMAX(residual, order=(p, d, q),
seasonal_order=(P, D, Q, 24)).fit(disp=-1)
#делаем предсказания линейной моделью
forecast = linear_model.predict(regression_feat[test_idx[:w_hours]])
start_ind = len(train_idx[lag_count * 24:])#индекс начала предсказания
# у нас y_true = y_pred - residual, поэтому вычитаем остатки
forecast -= arima_model.predict(start = start_ind, end = start_ind + 23)
actual = full_data[1231][test_idx[:w_hours]]
errors_rmse.append(rmse(actual.values, forecast))
errors_mape.append(mean_absolute_percentage_error(actual.values, forecast))
Посмотрим на результаты.
print(np.mean(errors_mape[2:]), np.mean(errors_rmse[2:]))
11.490083216107912 80.45325307865804
Видно что качество не стало лучше.
Было - Score: 11.880454673188066 80.4179657335404
Посмотрим на скор глазами.
errors_mape
[59.32325553829023, 12.294116607224375, 11.372000961378355, 9.84529860261997, 20.627729929000825, 10.05814308300142, 8.417277936190775, 10.128144896850317, 11.222562732876725, 10.249507586944913]
Видно, что в целом ошибка составляет примерно 10%, но например бывают дни, в которые модель ошибается на все 20%. (Не считая первые сплиты с маленьким трейном). Скорее всего, это является следствием того, что количество поездок на такси является очень шумным показателем, в том плане, что на него влияет очень много факторов, таких как, например, погода/пробки/праздники, которые не учитываются нашей моделью, но очень сильно могут влиять на спрос такси (значительно завышая/занижая его). Было бы интересно попробовать добавить подобные признаки, такие как предсказание синоптиков, в линейную модель).
Так что вполне возможно, что подобная точность в 11,5% является приемлимой.
Предсказание на один из последних дней выглядит как то так.
plt.plot(actual.values, label = 'y')
plt.plot(forecast.values, label = 'prediction')
plt.legend()
<matplotlib.legend.Legend at 0x1e2000839e8>
В данной работе мы решали задачу предсказания числа поездок такси в одном из центральных районов Нью-Йорка.
Выводы которые мы сделали в ходе знакомства с данными.
1. Основня часть поездок на такси в Нью-Йорке осуществляется из Манхэтенна.
2. Число поездок обладает, как дневной, так и недельной сезонностью. Это хорошо видно из графика.
Также, мы можем предположить, что есть и годовая сезонность.
Потому что число поездок в майе в среднем меньше (это мы проверили несколькими критериями), чем
в предыдущие более холодные месяцы. Это кажется вполне логичным, т.к чем теплее, тем более вероятно человек выберет
пройтись пешком, вместо того чтобы брать такси.
3. Поскольку ряд обладает сзонностью, то он не является стационарным. Однако, критерий Дики-Фулера, отвергает гипотезу о нестационарности.
4. Также, мы показали, что распределение числа поездок двумодальное, т.е люди катаются либо мало, либо много.
И оно, очевидно, не нормальное.
5. С помощью множественной проверки гипотез, мы показали что средние для различных месяцев отличаются.
Средние равны, только для 2 и 3, 2 и 4 месяца.
Для предсказания ряда, мы использовали линейную модель.
Мы сделали некоторые регрессионные признаки (день недели, час дня, лаги, гармонические признаки).
Проверяли качество на кросс-валидации.
В итоге удалось получить качество MAPE = 11,88%
Остатки от линейной модели оказались стационарными, несмещенными и похожими на белый шум.
Также мы попробовали обучить SARIMAX на остатках от линейной модели. Это не дало значительного прироста в качестве (удалось достигнуть MAPE = 11,49%).
Если смотреть на качество на разных фолдах (в ходе кросс-валидации), то можно увидеть, что есть дни в которые качество значительно проседает. Скорее всего это следствие того, что целевая переменная - число поездок на такси, довольно зашумлена, т.к она зависит от многих случайных факторов (погоды, пробки и т.д). Поэтому, возможно, что для первого приближения это неплохое решение.
Варианты улучшения:
1. Добавить, например, погодные признаки, индикаторы праздников и т.д
2. Можно рассмотреть в качестве признаков число поездок в соседних ячейках, т.к скорее всего они сильно коррелируют.
Это могло бы дать нам еще больше информации.
В данной задаче, решение с помощью линейной модели мне кажется оптимальным в плане простоты построения модели и настройки параметров. Ведь, изначально задача состояла в предсказании ряда для каждой из ячеек, и вариант, например, с построением модели ARIMA(или любой другой модели в которой надо кропотливо настраивать параметры) для каждой из 2500 ячеек кажется очень сомнительным.
P.S. За возможные грамматические ошибки прошу прощения.